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ABSTRACT 

We study a class of 2D non-equilibrium Ising models based on competing dynamics induced by 
contact with heat-baths at two different temperatures. We make a comparative study of the non- 
equilibrium versions of Metropolis, heat bath/Glauber and Swendsen-Wang dynamics and focus on 
their critical behavior in order to understand their universality classes. We present strong evidence 
that some of these dynamics have the same critical exponents and belong to the same universality 
class as the equilibrium 2D Ising model. We show that the bond version of the Swendsen-Wang 
update algorithm can be mapped into an equilibrium model at an effective temperature. 
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1.- Introduction 

Despite attempts at constructing a rigorous theory for non-equilibrium statistical mechanics, 
there is still no formalism to parallel the one which exists for equilibrium systems. As a result, 
there are few analytical methods with which to deal with non-equilibrium systems. In general, non- 
equilibrium systems display rich and complex behavior such as phase separation, pattern formation, 
turbulence M, 0], and it is therefore useful to first study simple model systems. 

There exist model non-equilibrium systems described by stationary distributions that are compar- 
atively easy to study. An open system maintained in a non-equilibrium steady state by an external 
temperature or density gradient is one example 0] . Another class of non-equilibrium steady states 
is obtained when the system is closed, and the dynamics is a local competition of two dynamics 
at different temperatures . It is this type of system that we address in this paper. Since these 
non-equilibrium systems display behavior qualitatively similar to equilibrium systems, such as phase 
transitions, it is important to ask if their critical properties and universality classes are the same. 
Flows of renormalized probability distributions and fixed points are normally independent of the 
details of a Hamiltonian. Perhaps these flows and fixed points are even independent of the existence 
of a Hamiltonian and a Boltzmannian distribution. If one were able to establish equivalences (or 
near equivalences) of universality classes between equilibrium and non-equilibrium models, then we 
will be able to answer many questions about the behavior of these systems without a complete 
formulation equivalent to the one for equilibrium systems. 

Simple non-equilibrium spin-flip stochastic systems have received considerable attention in the 
literature over the last 10 years. Reviews and general discussions about non-equilibrium phase 
transitions and stationary states can be found in refs. |Q, |5| g, 0|. Studies of driven diffusive systems 
can be found in Wang et al M , Marro et a/ S , Garrido et al U^ and Grinstein et al [0 . 

Grinstein et al [ p2[ studied the statistical mechanics of probabilistic cellular automata using time- 
dependent Ginzburg-Landau theory. They suggested that any non-equilibrium spin-flip dynamics 
with up-down symmetry belongs to the same universality class as the equilibrium Ising model. 
Their argument is based on the observation that under the renormalization group in rf = 4 — e, the 
dynamical fixed point of the Ising model is stable with respect to all additional analytic terms which 
preserve the lattice geometry and the spin up-down symmetry. 

Kanter and Fisher fl^ analyzed the existence of ordered phases in stochastic Ising systems with 
short-range interactions. They found that the existence of universality for those systems might de- 
pend on the details of the interaction and casted doubts on the general applicability of the argument 
of Grinstein et al |1^ . 

A two-temperature Glauber Ising model was introduced by Garrido, Labarta and Marro |14| (we 
will refer to this model as GLM) to investigate stationary non-equilibrium states. They obtained a 
mean field solution and performed some Monte Carlo simulations on 2D lattices. They found critical 
behavior qualitatively similar to the equilibrium case. 

The non-equilibrium behavior of competing dynamics such as spin-flip vs spin-exchange has been 
studied by Garrido et al [|l5| using hydrodynamic macroscopic equations and Monte Carlo data, and 
by Wang and Lebowitz p6[ using a Monte Carlo renormalization group method. They found evidence 
of equilibrium Ising behavior and Ising-like exponents. 

Tome et al n% studied the GLM model for the case when one of the temperatures is negative. 
They used a dynamical pair approximation to analyze antiferromagnetic steady states and obtained 



the corresponding phase diagram. Their mean field renormahzation group calculations show evidence 
in favor of equivalence with the equilibrium Ising universality class. 

Blote et al |1^ studied a model similar to GLM in which each sublattice is in equilibrium at a 
different temperature. They performed Monte Carlo simulations and found strong evidence that the 
model belongs to the equilibrium Ising universality class. 

Marques p9| , EOJ used a mean field renormahzation group calculation to obtain a phase diagram 
and calculated the exponent h' for the GLM model. Her results compared well with the equilibrium 
Ising values. Later she extended this technique to two different three state systems which retain 
the up-down symmetry [ pl| to test Grinstein et al. conjecture [n2|, and found good agreement with 
equilibrium Ising exponents. 

Recently de Oliveira p4 analyzed the isotropic majority-vote non-equilibrium model by Monte 
Carlo and finite size scaling. He found very good agreement for the critical exponents between 
this model and the equilibrium Ising model, and also for Binder's cumulant. In a separate paper 
de Oliveira et al p3| studied a family of non-equilibrium spin models with up-down symmetry 
parameterized by a Glauber-like transition rate. They also found good evidence that the critical 
exponents for this family (except for the limit case of the voter model) are the same as the equilibrium 
Ising model. 

Considerable numerical and analytic evidence has been accumulating in favor of universality and 
equilibrium Ising exponents for some of these models but the results are not definite yet. Most ana- 
lytical methods used, such as mean-field RG, are of an approximate nature. Monte Carlo simulations 
have focused mainly on qualitative behavior and the calculations of critical exponents have not been 
carried out with high resolution. To better understand and test this equivalence we have under- 
taken a detailed, high resolution Monte Carlo study of two-temperature Ising models and explore 
the question of universality using different local and non-local update dynamics. 

The paper is organized as follows. In Section 2 we will describe the different two-temperature 
non-equilibrium models that are the subject of this study. Section 3 presents detailed results for the 
Metropolis non-equilibrium dynamics including critical exponents and cumulant behavior. Section 
4 focuses on the Swendsen-Wang p4| non-equilibrium dynamics. Section 5 contains a comparative 
study of all the dynamics. Extension to many temperature model is discussed in Section 6 and the 
conclusions are presented in Section 7. 

2. Two-temperature non-equilibrium dynamics 

We begin with the two dimensional Ising model on the square lattice with Hamiltonian 

H = -JI3J2'^^'^J (1) 

where /3 is the inverse temperature and J the coupling. A dynamics for the model can be described 
in terms of a time dependent probability distribution P{a, t) which evolves according to a master 
equation. 



where W{a' *— a) is the transition rate from configuration cr to a'. We are interested in stationary 
probabiUty distributions, 

^ = 0, P(.,t).P(a) (3) 

In the case of equihbrium systems -P(cr) is not only stationary but has the form of a Boltzmann 
distribution parameterized by the inverse temperature /3, 

P{a) = Z-'^ exp(-/3iJ) (4) 

where Z ~ X^o- '2^P(~/^-^('^))- ^^^ our case we are interested in stationary non-equihbrium distri- 
butions produced by the local competition of equilibrium dynamics at different temperatures. The 
usual condition to obtain an equilibrium Monte Carlo dynamics is to make W{<7' ^- a) obey detailed 
balance (note that imposing detailed balance on the W^ is a sufficient but not necessary condition), 

W{a' ^a)P{a) ^W{a ^a')P{a'). (5) 

In the two-temperature model one considers a composite rate W{ij' ^- cr), 

W{a' <- a) =pWi{a' ^ a) + (1 -p)W^2(cr' <- a) (6) 

with competing Wi and W2 ■ At each time step the transition probability will be chosen at random 
to be Wi with probability p or W2 with probability {1 — p). W\ and W2 individually correspond to 
equilibrium transition rates obeying detailed balance with respect to temperatures f5\ and /32, «-e. 

W,{g ^ a') ^ eft^(") 

At each time step the dynamics obeys detailed balance locally with respect to /3i or P2 and the spins 
act as if in instantaneous contact with one of two heat baths. The overall effect is to produce a 
non-equilibrium dynamics which reduces to the equilibrium model when /?i = /?2- From the master 
equation one can prove that for a combined dynamics of this sort there is indeed a stationary regime 
given by the condition Q, 

p{aWi) + {l-p){<TW2)^0 (8) 

The induced global probability distribution P{<7) does not, in general, correspond to a local known 
Hamiltonian and it depends on the details of the dynamics, i.e. the particular choice of Wi and W2. 
This is in contrast to equilibrium simulations where the Boltzmann distribution is independent of 
the particular choice of W. The combination of Wi and W2 produces a stationary state analogous 
to a system being driven by an external potential. The ensemble of stationary configurations ex- 
hibits physical properties qualitatively similar to equilibrium (i.e. ordering, cluster formation, phase 
transitions, etc.). This is therefore one of the simplest ways to generate non-equilibrium models 
from equilibrium ones. In order to investigate the properties of stationary distribution on the dy- 
namics we study three different update algorithms: Metropolis, Glauber (or equivalently heat-bath) 
and Swendsen-Wang. Furthermore each of these dynamics has a "bond" and "spin" version as we 
explain below. 

We start by defining the GLM dynamics flj]. The transition rate for spin i takes the standard 
Glauber form, 

W^ = ^a[l-ta.nh{Jf3^a^ ^ a^)] (9) 

b-*l=i 



where j3i is chosen to be equal to /3i with probabihty p and j32 with probabihty (\ — p). In one 
dimension this dynamics is always equivalent to an equilibrium model at an effective temperature 
/3off given by Q, 

tanh(2J/3cff) = _ptanh(2J/?i) + (1 - p) tanh(2J/32) (10) 

Similarly "bond" dynamics is defined by 

t^z = ^a[l - tanh(jCT, ^ /3jaj)] (11) 

where (3j is (3i with probability p or (32 with probability (1 — p) and in this way the temperature is 
selected independently for each bond. These transition rates with a = 1 also define the heat-bath 
algorithm, i.e. the two algorithms are equivalent and we only need to discuss one. Henceforth we 
shall refer to this as the Glauber dynamics. 

Metropolis non-equilibrium dynamics are defined in a similar way. The relevant acceptance factor 
for spin i is given by, 

where, as in the Glauber case, /3 is either Pi or /?2. In the bond version we have, 

A.^e'^-'"'^^"'''' (13) 

where each bond j is chosen independently with temperature /3i or (32 ■ 

Finally we define non-equilibrium Swendsen-Wang |24[| spin and bond versions. In the bond 
version the percolation probability for bond ij is chosen to be, 

where /3y is /3i with probability p or /?2 with probability (1 — p). The subsequent percolation, 
cluster finding and flipping steps are done in the same way as in the original Swendsen-Wang 
dynamics. In section 4 we show that this dynamics can be mapped to an equilibrium system at an 
intermediate "effective" temperature /?off- The spin version is defined similarly except that the four 
bonds contributing to the update of each red (or black) site are chosen to be at the same temperature 
(3i or (32. 

All these algorithms can be implemented very efficiently on a parallel computer such as the CM - 
5. In this paper we present detailed results for the Metropolis-spin and Swendsen- Wang-bond cases 
and make some comparisons with the other dynamics. 



3. Metropolis non-equilibrium dynamics 

We performed a careful investigation of the spin version of the Metropolis dynamics. Our update 
algorithm is parallel, so we simultaneously update all the red/black sites on the (checkerboard) 
lattice. For each sublattice site i we choose temperature Pi or P2 independently using a uniformly 
distributed random number. Then we compute the change of energy with the spin fiipped. 



A^;, = 2a,^CTj. (15) 



If Aii"; is less or equal to zero then the flip is always accepted, otherwise it is accepted with probability 
expi-PAE,). 

To address the question of the existence of an equivalent equilibrium system we consider the 
local transition rate for the combined dynamics, 



e 



-fef A£ ^ p,-P.AE ^ (^ _ p)e-&AB^ (^g) 



This equation has to be satisfied for all values of AE, however, only the cases for AE > are 
temperature dependent^. In one-dimension there is only one relevant case {AE = 4) and the 
equation is satisfied with a (3cS equal to 

/?eff = -\ \og[pe-^^^ + (1 - p)e-^^^], (17) 

as was reported by Garrido et al [[l4| . In two-dimensions one has to satisfy two equations (for 
A£' = 4 and 8), 

pe-^'^i + (1 - p)e-^^^ = e-4^=//, (18) 

pe^^f'^ + {I - p)e-^f'' = e'^'^'ff. (19) 

For fixed p, /3i and P2 each equation has the solution 

Ptff = -lloglpe-^^^ + il-p^e-^P^] (20) 

Plff = -ilog[pe-8ft^(^_p)g-8fe] (21) 

which are shown in Fig. 1 as a function of /32 for fixed f3i — 0.35 and p = 0.5. One can see that only 
for Pi — (32 is /3^ff — Peff- This shows that for the Metropolis spin algorithm the two-temperature 
system cannot be described by a nearest neighbor effective Hamiltonian. The same is true for a 
triangular lattice with three (six) nearest-neighbors as there are two (three) independent equations. 

Monte Carlo simulations show that the configurations generated over time are qualitatively sim- 
ilar to equilibrium configurations and the system exhibits an Ising-like phase transition. We have 
explored the critical behavior of this dynamics on 16^, 32^, 64^ and 128^ lattices. To locate the 
critical point we set p — 0.5 and Pi — 0.35, and searched for the transition as a function of /32. Our 
best estimate for the critical point is P2 = 0.6372(5). The scaling region around this critical point 
appears to be rather narrow, and the calculation of critical exponent is therefore quite difficult. We 
choose this particular set of values for /3i and P2 as they are far from the Ising critical point {Pi — P2 
= 0.440678), and hope that any non-equihbrium effects will be manifest. We computed the critical 
exponents using finite size scaling and Binder's cumulant analysis. We accumulated measurements 
over 5 x 10*^ for 16^, 10^ steps for 32^ and 64^ lattices and 9 x 10^ for 128^ 

The data for magnetization and susceptibility for the four different lattice sizes are given in Table 
1 and plotted in Figs. 2 and 3. These figures show a change in the curvature between P2 = 0.6370 
and P2 = 0.6372. On the basis of these data we estimate that the critical coupling is P2 = 0.6371(1). 
The slopes give an estimate for the exponents P/v and 7/1^, assuming that the corrections to the 
leading finite size scaling forms m ~ L^^^/^ and x ^ L^^^ ^'^'^ negligible. The results are 

P/v = 0.122(2) (22) 

^/u = 1.73(2), (23) 



^We thank R. Swcndsen for bringing this to our attention 



where the errors are deterrained as foUows. We first compute the statistical error in the magnetization 
and susceptibihty for each data point. The error in the exponents is then obtained from the mean 
square fit to a straight fine on a fog- log plot. 
The data for Binder's cumulant pS, 

C/-l-^^^ (24) 

are also given in Table 1 and plotted in Figure 4 as a function of j32- The value of this cumulant at 
the critical point is conjectured to be a universal number independent of lattice size |2^. Figures 
5a and b show the cumulant values for the two pairs of lattice sizes, (32^, 64^) and (64^, 128^), at 
different temperatures around the critical point. The solid straight is the best fit to the data and 
the dashed line in the figures corresponds to C/l/2 = Ul ■ The point of crossing of these two lines 
gives an estimate of the critical point. The crossing takes place at 

C/*(32vs64) = 0.610(2) (25) 

C/*(64vsl28) = 0.605(10) (26) 

and the corresponding estimates for /Sj are 0.6370(2) and 0.6360(12) respectively. These values are 
consistent with the estimate from finite size scaling given above. Our longest runs were done at 
(32 = 0.6372, at which temperature our estimates 

[/i6 = 0.611(1) (27) 

C/32 = 0.611(2) (28) 

C/64 = 0.612(6) (29) 

C/128 = 0.615(7) (30) 

compare very well with the value U* ~ 0.611(1) computed by Bruce |E6| for the 2D equilibrium 
Ising model, and C/g4 = 0.611(5) independently by us. Thus, we shall henceforth call /32 = 0.6372 
the critical point. 

/,From the scaling of the cumulant one can obtain an estimate of ly, 



dUL 



dp 



= L^'^'GiL^'^e) (31) 

/3c 



where e = {(3c — (3)/ (3. To do this we first compute the slope dUL/d(3 for L = 32, 64 and 128 using 
the data points near the critical point and then fit these values versus L using the above expression 
to obtain an estimate for ly. We find, 

i/ = 0.99(5). (32) 

Another estimate can be computed from dllL/dllL' in the critical region because this quantity should 
scale as 

Using the more precise data shown in Fig. 5a for L = 32 and 64 we obtain v = 0.95(8). The 
agreement between these results and the values for the equilibrium Ising model is quite good (j/ = 1, 
(3/v = 0.125, 'y/v = 1.75, U* = 0.611(1)) and provides strong evidence for their equivalence. 



To further confirm this equivalence we measured the probabihty distribution for the magnetiza- 
tion and energy at /3i = 0.35 and (32 = 0.6372 and compare them with those for the equihbrium 
model in Figs. 7a and b. The agreement is somewhat better for the magnetization than for the 
energy and qualitatively the distribution functions are equilibrium-like. The data are slightly more 
disordered than those for the critical equilibrium model, suggesting that /3| may be slightly larger 
than 0.6372. 

The bond version of the Metropolis dynamics has not been studied as extensively, and we post- 
pone its discussion until Section 5 where we compare the different dynamics. 

4. Swendsen-Wang non-equilibrium dynamics 

We start with the bond contribution to a global probability distribution, 

Py- = pe-P^i^-'^^'^^) + (1 -p)e-'^^(i-''''"^) (34) 

If we assume that Pij is always equivalent to an equilibrium distribution with coupling /?off , i.e. 

P^^ = pe-/3i(i-'^.'^.) + (1 _ p)e-/32(i^<..<.,) ^ g-/3c«(i-o-,a,)^ (35^ 

then there exists a solution satisfying this equation for the two possible values of the bond energy, 
pe-2ft + (i_p)e-2fe ^ e-2/3eff^ [(t, = -a^] (36) 

p+{l-p) = 1, [(T^^CT.j\. (37) 



The solution is. 



-2/3oft _ g-2/3: 



P^— m I^mT' (38) 



e" 



-2/3: 



or equivalently, for a set of values /3i , /32 and p, there always exists a P^s given by, 

/3eff = -i log[p(e-2^i - e-^P-) + e-^/^^] (39) 

that corresponds to an Ising equilibrium system. This implies that this dynamics is nothing more 
than equilibrium dynamics in disguise. If we set /?cff = Pc we find lines of critical points in the 
/3i — /32 plane given by 

/?! = -ilog[l(e-2^= - (1 -p)e~^P^)]. (40) 

2 p 

A set of these critical lines is shown in Fig. 6 for different values of p. Our simulations confirm this 
equivalence. Notice that for some extreme values of /3i {(32) there is no positive value of (32 (/3i) that 
yields critical behavior. 

The reason for the equivalence is that for the Swendsen-Wang bond dynamics there are, in- 
dependent of the number of spatial dimensions, only two equations of constraint as each bond is 
independently in contact with the heat bath. One is the trivial condition p + (1 — p) = 1 and the 
second is the desired result given in Eq. p9. In a different context, a similar analysis of different ways 
to satisfy local equations for block percolation in equilibrium systems has been made in ref |g^ . 

Lastly, we have measured autocorrelation times at criticality with (3i = 0.1, (32 — 2.318, and 
p — 0.5 (/3off — 0.440687) to see if the stochastic choice of temperatures accelerates the decorrelation 



process. We find that the autocorrelation times are comparable to the values for the standard 
equilibrium Swendsen-Wang values. This is as expected since the two models are locally equivalent. 

5. Other dynamics and comparisons 

We have investigated the spin and bond version of Glauber dynamics for only three combinations 
of /3i and f32 using 32^ lattices. The ensemble size in these runs is ^ 50, 000 update sweeps, 
significantly smaller than for Metropolis or Swendsen-Wang dynamics. In Figs. 8a, 8b, and 8c 
we compare the equation of state {(E) vs {rn?)) for the various dynamics for the three different 
combinations of /3i and /32. Figure 8a shows the results for /3i — 0.4, (32 — 0.6, Figure 8b corresponds 
to /3i = 0.1, 132 = 2.318 and Figure 8c to /3i = 0.2, P2 = 0.738464. These three sets of temperatures 
were chosen such that the corresponding /3eff for Swendsen-Wang bond dynamics corresponds to a 
cold iPeff == 0.49) , critical {Peff = 0.4406868) and hot {Peff = 0.40) system respectively. All the 
simulations were done with p = 0.5. The dashed line corresponds to the result for the equilibrium 
Ising model obtained on L = 32^ lattices. The errors in the data are roughly equal to the size of the 
symbols. 

^Froni Figs. 8a-c it is clear that the equation of state depends very sensitively on the dynamics 
and there is a large spread in the results for all three choices of temperatures. In all cases, except for 
the Swendsen-Wang spin version, the results lie very close to the line characterizing the equilibrium 
Ising model. The deviations from the equilibrium results in all three versions of the 2-teinperature 
spin dynamics are in the direction of a more ordered system. The situation is reversed for the 
bond dynamics; compared to the equilibrium values the two-temperature results are less ordered. 
Qualitatively, the data with the various dynamics show the following ordering: Metropolis spin, 
Glauber spin, Metropolis bond, Glauber bond, and Swendsen-Wang bond. This pattern is also 
shown in Figs. 9a and b where we give the magnetization and energy probability distributions for 
/3i = 0.35 and (32 = 0.6372 (a critical point for the Metropolis spin dynamics). All cases appear to 
have the same functional form but their position and amplitudes are rescaled. The magnetization 
probability distributions are Gaussian near the peak but have long tails toward m = 0. Further 
work is needed to explore the possibility that all of these probability distributions are just rescaled 
forms of the equilibrium distribution and correspond to different /3eff's 

6. Three- Temperature Model 

To study the general case of the competition of many temperatures we extended the analy- 
sis to the three-temperature Metropolis spin dynamics. The three temperatures are chosen with 
probabilities pi, p2 and p^: 

W{a' ^ a)^piWi{a' ^ a) + P2W2W ^ a) + psWaia' ^ a). (41) 

We fixed /3i = 0.35, /^a = 0.6372, andpi = P2 — Ps = 1/3 and varied P2 about the equilibrium critical 
value P2 = 0.4406868. We expected the system to display critical behavior for P2 = 0.4406868. 

The system effectively displays critical behavior in that region as can be seen in Fig. 10 where 
results for Binder's cumulant are shown for (32 in the range [0.424, 0.456. Furthermore, the prob- 
ability distribution functions for (32 = 0.4406868 match the ones for the two-temperature model 



discussed in Section 4 and the equilibrium model at criticality as can be seen in Figs. 11a and b. 
We do find an apparent narrowing of the critical region compared with the two-temperature models 
and the statistical quality of the data are not accurate enough to measure the exponents. 

Based on the study of the 3 temperature model we make the following conjecture for the Metropo- 
lis spin dynamics. A model in contact with an arbitrary number of heat-baths will display Ising 
critical behavior provided each temperature or a pair of them is tuned to the critical value. As 
the number of pairs of temperatures increase the critical region becomes narrower, making the 
measurement of critical exponents and the study of critical behavior more difficult. 

7. Conclusions 

We present high statistics results showing that for the Metropolis spin dynamics the stationary 
states produced by the two-temperature model are very similar to equilibrium states. On basis of the 
agreement of the critical exponents and Binder's cumulant we conclude that the two-temperature 
Metropolis spin dynamics is in the same universality class as the Ising model. We also show that the 
bond version of the two-temperature Swendsen-Wang dynamics can be mapped into an equilibrium 
Ising model at an intermediate effective temperature. Thus, for these cases our results agree with the 
conjecture of Grinstein et al. that any non-equilibrium spin-flip dynamics which preserves up-down 
symmetry belongs to the same universality class as the equilibrium Ising model. Assuming that the 
system evolves into an equilibrium distribution after some thermalization steps, one can measure the 
critical exponents and flow of renormalized couplings using the Monte Carlo renormalization group 
method. We hope to investigate this possibility in the future. 

Our results for the Metropolis bond, Glauber spin and bond, and Swendsen-Wang spin dynamics 
are qualitative. Further work is required to conflrm that they too belong to the same universality 
class as the equilibrium Ising model. 

We have extended the two-temperature critical Metropolis spin dynamics to the three temper- 
ature case. We find that the system shows critical behavior when the third temperature is tuned 
to /3 = 0.4406868. Based on this we conjecture that the Ising critical behavior is preserved as long 
as one adds pairs of temperatures that are, by themselves, critical. The critical region appears 
to become narrower as the number of pairs of temperature values are increased and the statistical 
quality of the data deteriorates. This makes the calculation of critical exponents and the study of 
the models more difficult. 
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Table I: Data for Susceptibility, Magnetization and Binder's cuniulant for the 2 temperature 
Metropolis spin dynamics (/3i — 0.35, p — 0.5) 



A 


L 


X 


\rn\ 


U 


0.6372 


16 


9.11(5) 


0.7043(3) 


0.611(1) 


0.6300 
0.6350 
0.6370 
0.6372 
0.6375 
0.6400 
0.6475 
0.6500 
0.6550 
0.6600 


32 
32 
32 
32 
32 
32 
32 
32 
32 
32 


33.5(3) 
32.0(4) 
30.2(8) 
30.6(3) 
30.6(3) 
30.0(6) 
27.2(4) 
26.3(6) 
25.1(4) 
23.3(7) 


0.629(2) 

0.639(2) 

0.647(3) 

0.6464(7) 

0.6463(8) 

0.651(2) 

0.666(1) 

0.671(2) 

0.679(1) 

0.689(2) 


0.603(3) 

0.608(4) 

0.611(11) 

0.611(2) 

0.611(3) 

0.613(7) 

0.619(3) 

0.621(5) 

0.624(4) 

0.628(5) 


0.6300 
0.6350 
0.6370 
0.6372 
0.6375 
0.6400 
0.6475 
0.6500 
0.6550 
0.6600 


64 
64 
64 
64 
64 
64 
64 
64 
64 
64 


114(4) 

112(5) 

103(2) 

101(2) 

101(2) 

96(5) 

76(5) 

68(3) 

64(5) 

56(3) 


0.566(4) 
0.580(4) 
0.592(2) 
0.594(2) 
0.594(2) 
0.605(2) 
0.635(4) 
0.644(3) 
0.655(3) 
0.669(2) 


0.600(13) 

0.605(12) 

0.611(4) 

0.612(6) 

0.612(5) 

0.616(3) 

0.629(14) 

0.633(11) 

0.636(11) 

0.640(4) 


0.6350 
0.6370 
0.6372 
0.6375 
0.6400 
0.6475 
0.6500 
0.6550 
0.6600 


128 
128 
128 
128 
128 
128 
128 
128 
128 


409(51) 
357(22) 
315(12) 
328(18) 
257(20) 
197(19) 
148(10) 
133(15) 
91(10) 


0.52(1) 

0.538(6) 

0.550(2) 

0.550(3) 

0.573(5) 

0.611(3) 

0.626(3) 

0.644(3) 

0.664(2) 


0.597(40) 

0.607(22) 

0.615(7) 

0.615(11) 

0.627(21) 

0.639(9) 

0.645(10) 

0.649(9) 

0.654(7) 
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0,6 — /3i = 0,35, p = 0,5 
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Fig. 1.- Plot of Pjjj and /3fj^ (see eq. |^ versus (32 (at /3i 
non-equilibrium Metropolis spin dynamics. 
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Fig. 2.- Scaling behavior of Magnetization as a function of /32 (with /3i = 0.35 and p = 0.5 held 
fixed) on different lattice sizes for the non-equilibrium Metropolis spin dynamics. 
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Fig. 3.- Scaling behavior of susceptibility as a function of f32 (with /3i — 0.35 and p = 0.5 held 
fixed) on different lattice sizes for the non-equilibrium Metropolis spin dynamics. 
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Fig. 4.- Cumulant values as a function of f32 (at /3i = 0.35 and p = 0.5) and on different lattice 
sizes for the non-equilibrium Metropolis spin dynamics. 
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Fig. 5a.- Cumulant values for 32^ versus 64^ lattices at temperatures around the critical point for 
the non-equilibrium Metropolis spin dynamics. 
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Fig. 5b.- Cumulant values for 64^ versus 128^ lattices at temperatures around the critical point 
for the non-equilibrium Metropolis spin dynamics. 
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Fig. 6.- Critical lines (see eq. koh for the Swendsen-Wang bond dynamics. 
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Fig. 7a.- Magnetization probability distributions P(|?7i|) for the Metropolis spin dynamics at 
/3i — 0.35, /32 = 0.6372, p = 0.5 and L = 32 (critical point). The probability distributions for the 
2D and 31) equilibrium Ising model are also shown for comparison. 
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Fig. 7b.- Energy probability distributions P{—E) for the Metropolis spin dynamics at /3i — 0.35, 
(32 = 0.6372, p = 0.5 and L — 32 (critical point). The probability distributions for the 2D and 31? 
equilibrium Ising model at criticality are also shown for comparison. 
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Fig. 8a.- Equation of state data for different dynamics with /3i = 0.4, (32 = 0.6, p = 0.5 and L = 32 
(^SM^.o„. = 0.49) 
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Fig. 8b.- Equation of state data for different dynamics with /3i = O.l, (32 = 2.318, p = 0.5 and 
L = 32 (/3fff '"'"'' = 0.4406868) 
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Fig. 8c.- Equation of state data for different dynamics with f3i = 0.2, (32 = 0.738464, p = 0.5 and 
i = 32 {(i^^""^" = 0.4) 
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Fig. 9a.- Magnetization probability distributions P(|?7i|) for the spin and bond versions of Metropo- 
lis, Glauber and Swendsen- Wang dynamics at (3i = 0.35, (32 = 0.6372, p — 0.5 and L — 32 (critical 
point). 
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Fig. 9b.- Energy probability distributions P(—E) for tlie spin and bond versions of Metropolis, 
Glauber and Swendsen- Wang dynamics at (3i — 0.35, (32 — 0.6372, p — 0.5 and i = 32 (critical 
point). 
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Fig. 10.- Cumulant values as a function of /?2 (at /3i = 0.35 and fS^ 
lattice sizes for the three-temperature Metropolis spin dynamics. 
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Fig. 11a.- Magnetization probability distributions P(|?7i|) for the two temperature and three tem- 
perature Metropohs spin dynamics at the critical point (L = 32). 



29 



0.02 



H 



^ 



0,015 



0.01 



0.005 



-0 



2 — temp Metro spin 
□ 3 — temp Metro spin 




1.2 



1.4 
- E 



1.6 



1.8 



Fig. lib.- Energy probability distributions P{—E) for tlic two temperature and three temperature 
Metropolis spin dynamics at the critical point {L = 32). 
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